% DHODAPsi2_fix.m - 修复数组越界问题的代码

% 根据工作区变量的信息，我们知道：
% I1, I2, I12 是 98x98 double
% k, l 是 1x98 double
% cols = 100, rows 可能也是100
% 但实际上我们只能使用到索引98

function [Psi2_val] = DHODAPsi2(signal)
    % 获取信号尺寸
    [rows, cols] = size(signal);
    
    % 初始化梯度计算
    I1 = zeros(rows-2, cols-2);
    I2 = zeros(rows-2, cols-2);
    
    % 计算梯度
    for k = 2:rows-1
        for l = 2:cols-1
            I1(k-1, l-1) = (signal(k+1, l) - signal(k-1, l)) / 2;
            I2(k-1, l-1) = (signal(k, l+1) - signal(k, l-1)) / 2;
        end
    end
    
    % 计算Matrix_flesh
    Matrix_flesh = max(abs(I1), abs(I2));
    
    % 初始化I2_peel
    I2_peel = zeros(rows-2, cols-2);
    
    % 设置I2_peel的边界
    I2_peel(:,1) = Matrix_flesh(:,1);
    I2_peel(:,end) = Matrix_flesh(:,end);
    I2_peel(1,:) = Matrix_flesh(1,:);
    I2_peel(end,:) = Matrix_flesh(end,:);
    
    % 设置内部值
    I2_peel(2:end-1, 2:end-1) = Matrix_flesh(2:end-1, 2:end-1);
    
    % 计算I12
    I12 = zeros(rows-4, cols-4);
    for k = 2:rows-3
        for l = 2:cols-3
            I12(k-1, l-1) = (I2_peel(k+1, l) - I2_peel(k-1, l)) / 2;
        end
    end
    
    % 计算Psi2_val，确保索引不会超出边界
    % 原始代码中的问题是索引超出了边界，我们需要调整计算范围
    Psi2_val = zeros(rows-4, cols-4);
    
    for k = 2:rows-3
        for l = 2:cols-3
            % 确保所有引用的索引都在有效范围内
            Psi2_val(k-1, l-1) = 2*signal(k, l)^2 ...
                - signal(k-1, l) * signal(k+1, l) ...
                - signal(k, l+1) * signal(k, l-1) ...
                + (2 * I1(k-1, l-1) * I2(k-1, l-1) - signal(k, l) * I12(k-1, l-1));
        end
    end
end 